Week 02 - Aarhus Environs: Introduction to Reading and Plotting Spatial Data

Introduction into GIS in R

There are a few key spatial packages available for Spatial Analysis in R, which you need to install as you progress through these exercises. The most basic are

  • sf for working with vector data
  • raster and terra for working with raster data

Task 1: Reading vector data

The sf package, created by Edzer Pebesma and colleagues, has dramatically simplified reading vector spatial data into R.

In this exercise you will read in three shapefiles (parks, playgrounds, and forests one point file and two polygon files) and one geojson (shelters) using st_read(). If you’ve read in the files correctly, you will see a standard R data frame except it will show some header metadata about the file and you’ll see a special geometry column which we will discuss later.

Instructions

  • Load the sf package.
  • All your datasets reside in the ‘data’ folder.
  • Import the forests shapefile (“forests.shp”).
  • Import the playgrounds shapefile (“playgrounds4326.shp”).
  • Import the parks shapefile (“parks.shp”).
  • Import the shelters geojson (“shelters.json”).
  • Use the head() function and identify the first few features of each layer.
# Load the sf package
___(___)

# Read in the forests shapefile
forests <- ___("../data/forests.shp")

# Read in the parks shapefile
parks <- ___(___)

# Read in the playgrounds shapefile
playgrounds <- read_sf("../data/playgrounds4326.shp", crs = 4326)

# Read in the shelters json
shelters <- ___(___)

# View the first few features of all layers
___(___)

Questions:

  1. How many features does each layer contain and what kind of geometry are they?
  2. What is the CRS value in these objects?

Solution

Reading layer `parks' from data source 
  `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\parks.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 47 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 561480.2 ymin: 6210157 xmax: 581216 ymax: 6237475
Projected CRS: ETRS89 / UTM zone 32N
Reading layer `forests' from data source 
  `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\forests.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 32 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 566833 ymin: 6212338 xmax: 582568.3 ymax: 6236083
Projected CRS: ETRS89 / UTM zone 32N
Reading layer `playgrounds4326' from data source 
  `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\playgrounds4326.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 51 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 9.990331 ymin: 56.03292 xmax: 10.30607 ymax: 56.26783
CRS:           NA
Reading layer `shelters' from data source 
  `C:\Users\Adela\Documents\RStudio\1_Teaching\cds-spatial\data\shelters.json' 
  using driver `GeoJSON'
Simple feature collection with 19 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 563459.7 ymin: 6212589 xmax: 578934.2 ymax: 6232898
Projected CRS: ETRS89 / UTM zone 32N
Simple feature collection with 6 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 570961.1 ymin: 6216613 xmax: 577431.9 ymax: 6232898
Projected CRS: ETRS89 / UTM zone 32N
  bookbar     type                               navn skov
1      Ja Shelters           Shelter i Gjellerup Skov <NA>
2      Ja Shelters     Shelter i H\xf8rhaven i skoven <NA>
3      Ja Shelters     Shelter i Lisbjerg gammel skov <NA>
4      Ja Shelters     Shelter ved Moesg\xe5rd Strand <NA>
5      Ja Shelters            Shelter i Mollerup Skov <NA>
6      Ja Shelters Shelter i H\xf8rhaven p\xe5 bakken <NA>
                                          mi_style mi_prinx
1 Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        1
2 Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        2
3 Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        3
4 Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        4
5 Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        5
6 Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        6
                  geometry
1 POINT (570961.1 6223407)
2 POINT (576307.6 6219521)
3 POINT (572490.7 6232898)
4 POINT (577431.9 6216613)
5 POINT (574609.9 6229691)
6 POINT (576455.1 6219324)
Simple feature collection with 6 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 566833 ymin: 6223092 xmax: 576573 ymax: 6232903
Projected CRS: ETRS89 / UTM zone 32N
  bookbar        type           navn skovnr areal_m2   min_x  min_y   max_x
1      Ja Skove_store    Åbyhøj Skov    290    32656 -224407 191537 -224289
2      Ja Skove_store  Brabrand Skov    240    39231 -225585 191149 -225221
3      Ja Skove_store    Årslev Skov    250    76629 -228610 191124 -228087
4      Ja Skove_store Gjellerup Skov    190   216967 -224931 190847 -224289
   max_y oprettet_a              oprettet_d rettet_af              rettet_dat
1 192017    az25000 2019-05-07 12:02:50.657   az25000 2019-05-07 12:02:50.657
2 191359    az25000 2019-05-07 12:02:50.657   az25000 2019-05-07 12:02:50.657
3 191374    az25000 2019-05-07 12:02:50.657   az25000 2019-05-07 12:02:50.657
4 191485    az25000 2019-05-07 12:02:50.657   az25000 2019-05-07 12:02:50.657
                                      mi_style mi_prinx
1 Pen (2, 1, 32768) Brush (2, 32768, 16777215)        1
2 Pen (2, 1, 32768) Brush (2, 32768, 16777215)        2
3 Pen (2, 1, 32768) Brush (2, 32768, 16777215)        3
4 Pen (2, 1, 32768) Brush (2, 32768, 16777215)        4
                        geometry
1 MULTIPOLYGON (((571023.4 62...
2 MULTIPOLYGON (((570135.1 62...
3 MULTIPOLYGON (((567335.9 62...
4 MULTIPOLYGON (((571154.8 62...
 [ reached 'max' / getOption("max.print") -- omitted 2 rows ]
Simple feature collection with 6 features and 10 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 567350.2 ymin: 6210157 xmax: 581216 ymax: 6237475
Projected CRS: ETRS89 / UTM zone 32N
  bookbar   type              navn areal_m2 oprettet_a              oprettet_d
1      Ja Parker Aarslev Møllepark    40022    az25000 2019-05-07 12:02:19.407
2      Ja Parker           Bananen     5667    az25000 2019-05-07 12:02:19.407
3      Ja Parker       Bavneparken    35611    az25000 2019-05-07 12:02:19.407
4      Ja Parker     Bækvejsparken    19846    az25000 2019-05-07 12:02:19.407
5      Ja Parker Engdalgårdsparken    81476    az25000 2019-05-07 12:02:19.407
6      Ja Parker       Forteparken    65751    az25000 2019-05-07 12:02:19.407
  rettet_af              rettet_dat
1   az25000 2019-05-07 12:02:19.407
2   az25000 2019-05-07 12:02:19.407
3   az25000 2019-05-07 12:02:19.407
4   az25000 2019-05-07 12:02:19.407
5   az25000 2019-05-07 12:02:19.407
6   az25000 2019-05-07 12:02:19.407
                                           mi_style mi_prinx
1 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215)        1
2 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215)        2
3 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215)        3
4 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215)        4
5 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215)        5
6 Pen (3, 1, 3825664) Brush (2, 11579392, 16777215)        6
                        geometry
1 MULTIPOLYGON (((567616.2 62...
2 MULTIPOLYGON (((581047.6 62...
3 MULTIPOLYGON (((567728.9 62...
4 MULTIPOLYGON (((574082.9 62...
5 MULTIPOLYGON (((575090.4 62...
6 MULTIPOLYGON (((575844.2 62...
Simple feature collection with 19 features and 6 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 563459.7 ymin: 6212589 xmax: 578934.2 ymax: 6232898
Projected CRS: ETRS89 / UTM zone 32N
First 10 features:
   bookbar     type                               navn skov
1       Ja Shelters           Shelter i Gjellerup Skov <NA>
2       Ja Shelters     Shelter i H\xf8rhaven i skoven <NA>
3       Ja Shelters     Shelter i Lisbjerg gammel skov <NA>
4       Ja Shelters     Shelter ved Moesg\xe5rd Strand <NA>
5       Ja Shelters            Shelter i Mollerup Skov <NA>
6       Ja Shelters Shelter i H\xf8rhaven p\xe5 bakken <NA>
7       Ja Shelters            Shelter p\xe5 Vestereng <NA>
8       Ja Shelters         Shelter i Moesg\xe5rd Skov <NA>
9       Ja Shelters           Dagshelter i H\xf8rhaven <NA>
10      Ja Shelters          Shelter i Brendstrup Skov <NA>
                                           mi_style mi_prinx
1  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        1
2  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        2
3  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        3
4  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        4
5  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        5
6  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        6
7  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        7
8  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        8
9  Symbol (34,8388608,14,"MapInfo Real Estate",1,0)        9
10 Symbol (34,8388608,14,"MapInfo Real Estate",1,0)       10
                   geometry
1  POINT (570961.1 6223407)
2  POINT (576307.6 6219521)
3  POINT (572490.7 6232898)
4  POINT (577431.9 6216613)
5  POINT (574609.9 6229691)
6  POINT (576455.1 6219324)
7  POINT (573438.3 6227501)
8  POINT (576374.9 6217600)
9  POINT (576498.6 6219401)
10 POINT (572113.1 6227183)
Simple feature collection with 51 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 9.990331 ymin: 56.03292 xmax: 10.30607 ymax: 56.26783
CRS:           NA
First 10 features:
   ogr_fid nr         lokatin                               adresse
1        1  1 A. H. Wingesvej    A. H. Wingesvej 13, 8000 Aarhus C.
2        2  2        Agernvej               Agernvej 73, 8330 Beder
3        3  3        Baunevej          Baunevej 15, 8361 Hasselager
4        4  4        Bellevue            Fortevej 109, 8240 Risskov
5        5  5    Beringsminde         Tranebærvej 44, 8220 Brabrand
6        6  6   Botanisk Have v. Peter Holms vej 17, 8000 Aarhus C.
7        7 61   Botanisk Have        v. Mønsgade 12, 8000 Aarhus C.
8        8  7        Byvangen         Grundvigsvej 55, 8260 Viby J.
9        9  8   Bækvejsparken               Bækvej 31, 8340 Malling
10      10  9       Egå Engsø                 Viengevej 7, 8250 Egå
                 X_tmstm                  geometry
1  2020-02-14 00:00:00.0 POINT (10.20414 56.17545)
2  2020-02-14 00:00:00.0 POINT (10.20171 56.06382)
3  2020-02-14 00:00:00.0 POINT (10.09647 56.10748)
4  2020-02-14 00:00:00.0 POINT (10.24744 56.19133)
5  2020-02-14 00:00:00.0 POINT (10.12097 56.16046)
6  2020-02-14 00:00:00.0 POINT (10.19203 56.16155)
7  2020-02-14 00:00:00.0 POINT (10.19378 56.15858)
8  2020-02-14 00:00:00.0 POINT (10.17067 56.13048)
9  2020-02-14 00:00:00.0 POINT (10.19102 56.03292)
10 2020-02-14 00:00:00.0 POINT (10.22041 56.21385)
                   name                                  long_name write  copy
PCIDSK           PCIDSK                       PCIDSK Database File  TRUE FALSE
netCDF           netCDF                 Network Common Data Format  TRUE  TRUE
PDS4               PDS4               NASA Planetary Data System 4  TRUE  TRUE
VICAR             VICAR                            MIPL VICAR file  TRUE  TRUE
JP2OpenJPEG JP2OpenJPEG JPEG-2000 driver based on OpenJPEG library FALSE  TRUE
PDF                 PDF                             Geospatial PDF  TRUE  TRUE
MBTiles         MBTiles                                    MBTiles  TRUE  TRUE
BAG                 BAG                 Bathymetry Attributed Grid  TRUE  TRUE
EEDA               EEDA                      Earth Engine Data API FALSE FALSE
OGCAPI           OGCAPI                                     OGCAPI FALSE FALSE
            is_raster is_vector   vsi
PCIDSK           TRUE      TRUE  TRUE
netCDF           TRUE      TRUE FALSE
PDS4             TRUE      TRUE  TRUE
VICAR            TRUE      TRUE  TRUE
JP2OpenJPEG      TRUE      TRUE  TRUE
PDF              TRUE      TRUE  TRUE
MBTiles          TRUE      TRUE  TRUE
BAG              TRUE      TRUE  TRUE
EEDA            FALSE      TRUE FALSE
OGCAPI           TRUE      TRUE  TRUE
 [ reached 'max' / getOption("max.print") -- omitted 62 rows ]
[1] "C:/Users/Adela/Documents/RStudio/1_Teaching/cds-spatial/Week02"

Well done, now you should see how easy it can be to read in shapefiles and you got your first taste of what an sf object looks like.


Task 2: sf objects are data frames

Spatial objects in sf are just data frames with some special properties. This means that packages like dplyr can be used to manipulate sf objects. In this exercise, you will use the dplyr functions select() to select or drop variables, filter() to filter the data and mutate() to add or alter columns.

Instructions

  • Load the dplyr and sf packages.
  • Use the filter() function from dplyr on the parks object to create a new data frame limited to parks greater than 5ha.
  • Use the nrow() function on your new object to determine how many parks greater than 5ha are in the dataset.
  • Use the select() function from dplyr to limit the variables in your over5ha dataset to just navn and areal_m2 and create a new data frame.
  • Use the head() function to check which variables exist in your new data frame. Does the data frame only have the navn and areal_m2 columns (the answer is no, why)?
# Load the sf package
___

# ... and the dplyr package
___

# Use filter() to limit to over5ha parks
___

# Count the number of rows
___(over5ha)

# Limit to navn and areal_m2 variables
over5_lim <- over5ha %>% ___(navn, areal_m2) 

# Use head() to look at the first few records
___(over5_lim)

Solution

[1] 17

Great! You can see why the sf package is so nice – your spatial objects are data frames that you can smoothly manipulate with dplyr. The number of parks over 5ha is 17 You may have noticed that when you used select the default is to keep the geometry column even if you didn’t explicitly list it as a column in select.

Task 3: Geometry is stored in list-columns

A major innovation in sf is that spatial objects are data frames. This is possible thanks, in part, to the list-column.

A list-column behaves, to a certain extent, like any other R column. The main difference is that instead of a standard value such as a single number, character or boolean value, each observation value in that column is a piece of an R list and this list can be as complex as needed. The list column allows you to store far more information in a single variable and sf takes advantage of this by storing all geographic information for each feature in the list.

In this exercise, you will convert the data frame to what’s called a tibble with tibble::as_tibble() (Note that dplyr::tbl_df() is now deprecated).

Instructions

  • Load tidyverse in your workspace.
  • Create a simple data frame df that includes a single column a using data.frame().
  • Add a list-column b to your data frame with the list() function.
  • Use head() to look at df.
  • Use as_tibble() to convert the data frame to a tibble and print it to the console. This is just for cleaner printing.
  • Pull out the third observation from columns a and b using base R (you’ll need square brackets like [3]).
# Create a standard, non-spatial data frame with one column
df <- ___(a = 1:3)

# Add a list column to your data frame
df$b <- ___(1:4, 1:5, 1:10)

# Look at your data frame with head
___(df)

# Convert your data frame to a tibble and print on console
___(df)

# Pull out the third observation from both columns individually
df$___[___]
df$___[___]

Solution

  a                             b
1 1                    1, 2, 3, 4
2 2                 1, 2, 3, 4, 5
3 3 1, 2, 3, 4, 5, 6, 7, 8, 9, 10
# A tibble: 3 × 2
      a b         
  <int> <list>    
1     1 <int [4]> 
2     2 <int [5]> 
3     3 <int [10]>
[1] 3
[[1]]
 [1]  1  2  3  4  5  6  7  8  9 10

You now have a better sense of what a list column is. You can see how it can be used to store far more information in a single variable than other types of columns. These list-columns are how sf stores detailed geographic information on each feature in a single record. Converting the data frame to a tibble is not necessary but a tibble can provide a better print out of the object.

Task 4: Extract geometric information from your vector layers

Sometimes you spatial data will come with attributes such as area, length, etc. Unless you have just created these yourself, you cannot be sure of their provenance and accuracy. They may have also changed in the course of your spatial transformations. If you want to know what is the actual area of your polygons, then you can extract it directly from the geometries.

There are several functions in sf that allow you to access geometric information like area from your vector features. For example, the functions st_area() and st_length() return the area and length of your features, respectively.

Note that the result of functions like st_area() and st_length() will not be a traditional vector. Instead the result has a class of units which means the vector result is accompanied by metadata describing the object’s units. As a result, code like this won’t quite work:

# This will not work
result <- st_area(parks)
result > 100000

Instead you need to either remove the units with unclass():

# This will work
val <- 100000
which(unclass(result) > 100000)

or you need to convert val’s class to units, for example:

# This will work
units(val) <- units(result)
which(result > val)
length(which(result > val))

Instructions

  • Check that your sf library is still active and the forests shapefile in memory
  • Compute the area of the forest units.
  • Create a histogram of the areas using hist() to quickly visualize the data spread.
  • Filter the forests object with filter() and limit to forests with unclass(areas) > 200000 (areas greater than 20 hectares).
  • Can you plot the geometry and their names?
  • Plot the geometry of the result with plot() and st_geometry().
# Compute the areas of the forests
areas <- ___(forests)

# Create a quick histogram of the areas using hist
___(___)

# Filter to forests greater than 200000 (square meters)
big_forests <- ___ %>% ___(___(___) > 200000)

# Can you plot the big_forests with their names?(hint: check the plotting section in Task 5 below)


# Plot just the geometry of big_forests
___(___(big_forests))

Solution

 [1] "Gjellerup Skov"               "Mollerup Skov"               
 [3] "Vestereng"                    "Tranbjerg Skov"              
 [5] "Vilhelmsborg Skov"            "Skødstrup Skov"              
 [7] "Hørret Skov"                  "Brendstrup Skov"             
 [9] "Viby Høskov"                  "Fløjstrup Skov"              
[11] "Kirkeskov"                    "Lisbjerg Skov, den gamle del"
[13] "Lisbjerg Skov, den nye del"   "Thorsskov"                   
[15] "Hestehaven"                   "Skåde Skov"                  
[17] "Havreballe Skov"              "Riis Skov"                   
[19] "Moesgård Skov"                "Mariendal Skov"

Excellent! Computing geographic information for your vector layers can be done with functions like st_area() and st_length(). As you saw in this exercise, these functions produce a result that can be used in additional calculations but you need to be careful because the result is a units object that requires a little additional processing like using unclass().

Task 5: Plot vector spatial objects

The function for making a quick map/plot is a function you are already familiar with, plot(). You can, for example, type plot(my_data) to see your spatial object. The default, though, may not be what you want. The plot() function, when applied to sf objects, will create a set of maps, one for each attribute in your data. Instead, if you want to create a map of a single attribute you can extract that attribute using, as an example, plot(my_data["my_variable"]).

Frequently you just want to plot the raw geometry with no attribute color-coding (e.g., adding county boundaries to a map of points). For this, you can use the st_geometry() function to extract the geometry and plot the result. You can either create a new object or you can nest st_geometry() within the plot() function.

Often, you also want to plot multiple spatial objects together to see if and how they relate. To do this, use the plot()function twice in sequence, separated by ; with the second instance containing the add=TRUE argument plot( ,add=TRUE). You can use the col argument to differentiate each layer by color.

Instructions

  • Use plot() to plot the forest data using all defaults.
  • Plot just the areal_m2 attribute of forests.
  • Create a new object that is just the geometry of the forests object with st_geometry().
  • Plot the geometry of the forests (the object you just created).
  • Plot the geometry of the forests and the parks together, making the parks pink and the forests green.

Solution

Well done! Yes, these plots are not super pretty but you can’t beat plot() for a quick look using few keystrokes. And remember you can use plot(st_geometry(geo_object)) to plot just the geometry of your object.

Task 6: Reading in raster data

The term “raster” refers to gridded data that can include satellite imagery, aerial photographs (like ortophotos) and other types. In R, raster data can be handled using the raster package created by Robert J. Hijmans or by the newly developed terra package by the same author. See which one you can get to work in your workspace.

When working with raster data, one of the most important things to keep in mind is that the raw data can be what is known as “single-band” or “multi-band” and these are handled a little differently in R. Single-band rasters are the simplest, these have a single layer of raster values – a classic example would be an elevation raster where each cell value represents the elevation at that location.

Multi-band rasters will have more than one layer. An example is a color aerial photo in which there would be one band each representing red, green or blue light (RGB).

Instructions

  • Load the raster or terra packages.
  • Load the elevation grid DNK_msk_alt.grd with the raster function and assign to elevation object
  • Read in the orthophoto image with terra:rast() (this is multi-band raster called “Aarhus_1m.TIF”).
  • Use the class() function to determine the class of each raster object you read in.
  • Use the nlayers() and nlyr() function to determine how many bands/layers are in each object.
  • Use the res() function to learn the resolution of the image
# Load the raster package
___(___)

# Read in the mound elevation single-band raster
elevation <- ____(___)

# Read in the orthophoto image multi-band raster
aarhus <- brick(_____)

# Get the class for the new objects
class(___)
class(___)

# Identify how many layers each object has
nlayers(___)
nlayers(___)

# Identify the resolution of each raster
___(elevation)
___(aarhus)

Questions:

  1. What are the dimensions (number of columns and rows) of the elevation raster?
  2. What is the resolution of aarhus orthophoto? How many layers does it contain and what do they represent?

Now you’ve learned how to read in single and multi-band rasters. You should have noticed, based on the nlayers() function, that the elevation object has a single layer and the aarhus object has five (well, four useful ones, if you plot it).

Solution

[1] "RasterLayer"
attr(,"package")
[1] "raster"
[1] "SpatRaster"
attr(,"package")
[1] "terra"
[1] 1
[1] 5
[1] 0.008333333 0.008333333
[1] 1 1

Task 7: Learn about your raster objects

Instead of storing raster objects in data frames, the raster\terra package stores spatial data in specially designed R classes that contain slots where the data and metadata are stored. The data and metadata can be accessed using a suite of functions. For example, the spatial extent (the bounding box) of the object can be accessed with extent()\ext(), the coordinate reference system can be accessed with crs() and the number of grid cells can be determined with ncell().

Instructions

  • You should have raster and terra packages loaded, and have the elevation layer and ortophoto image layer for Aarhus in memory.
  • Use the extent() and ext() functions to get the extent of the elevation and aarhus layer.
  • Use the crs() function to get the coordinate reference system of aarhus and elevation.
  • Use the ncell() function to determine how many grid cells are in the elevation layer and the aarhus layer.
# Get the extent of the elevation and aarhus object
___(___)

# Get the CRS of the aarhus and elevation object
___(___)

# Determine the number of grid cells in both raster objects
___(aarhus)
___(elevation)

Solution

class      : Extent 
xmin       : 8 
xmax       : 15.3 
ymin       : 54.5 
ymax       : 57.8 
SpatExtent : 573334, 576668, 6223334, 6226668 (xmin, xmax, ymin, ymax)
[1] "PROJCRS[\"ETRS89 / UTM zone 32N\",\n    BASEGEOGCRS[\"ETRS89\",\n        DATUM[\"IRENET95\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1,\n                    ID[\"EPSG\",9001]]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]],\n    CONVERSION[\"Transverse Mercator\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=WGS84 +no_defs 
WKT2 2019 representation:
GEOGCRS["unknown",
    DATUM["Unknown based on WGS 84 ellipsoid",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1],
            ID["EPSG",7030]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]] 
[1] 11115556
[1] 346896

Great work! Although rasters are not stored as data frames, the metadata can easily be extracted using functions like extent() or ext(), crs() and ncell().

Question 3: It makes sense that the extents and number of cells for Danish national elevation data and Aarhus city would be different because their sizes diverge, but why are the units so different? What are these units?

Task 8: Plot your raster object

Similar to what you saw in the exercises related to vector objects it’s often useful to quickly look at a map of your raster objects with the plot() function.

The raster package has added useful methods for plotting both single and multi-band rasters. For single-band rasters or for a map of each layer in a multi-band raster you can simply use plot(). If you have a multi-band raster with layers for red, green and blue light you can use the plotRGB() function to plot the raster layers together as a single image.

Instructions

  • Plot the elevation raster with the plot() function, it is a single-band raster.
  • Plot the aarhus object with the plot() function, it is a multi-band raster.
  • Plot the aarhus raster with plotRGB() to see all layers plotted together as a single image. Check out the optional stretch argument in this function, and set it to linear to brighten the image up a little
  • Try to plot the two images together (on top of one another). Can you do it?
# Plot the elevation raster (single raster)
___

# Plot the aarhus raster (as a single image for each layer)
___

# Plot the aarhus raster as an image
____

# Plot the two images together (on top of one another). Hint: it only works if their CRSes are compatible. Return to this after Task 9 and 10 if the CRSes are different.
____

Solution

Nice work! As you can see, the plot() function can be used to plot single layers while the plotRGB() function can be used to combine layers into a single image. Plotting two raster objects of different extents, resolution and CRS can not be done without additional wrangling

Task 9: Vector and raster coordinate systems

In order to perform any spatial analysis with more than one layer, your layers should share the same coordinate reference system (CRS) and the first step is determining what coordinate reference system your data has. To do this you can make use of the sf function st_crs() and the raster\terra function crs().

When the geographic data you read in with sf already has a CRS defined both sf and raster\terra will recognize and retain it. When the CRS is not defined you will need to define it yourself using either the EPSG number or the proj4string.

Instructions

  • Ensure the packages sf and raster\terra and the objects forests and playgrounds and aarhus and elevation are loaded in your workspace.
  • Use st_crs() to identify if a CRS exists and what it is for the playgrounds and forests objects.
  • Use the st_crs() function to define/assign a CRS to the playgrounds object, utilising the EPSG number 4326 as that is the original projection (now lost).
  • Use crs() to identify if a CRS exists and what it is for the aarhus and elevation objects.
  • Do not use the crs() function to define/assign a CRS to the elevation object as it already has a CRS of its own and renaming it won’t change the properties of the file.

Solution

# Determine the CRS for the elevation and playgrounds vector objects
st_crs(forests)
Coordinate Reference System:
  User input: ETRS89 / UTM zone 32N 
  wkt:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        ENSEMBLE["European Terrestrial Reference System 1989 ensemble",
            MEMBER["European Terrestrial Reference Frame 1989"],
            MEMBER["European Terrestrial Reference Frame 1990"],
            MEMBER["European Terrestrial Reference Frame 1991"],
            MEMBER["European Terrestrial Reference Frame 1992"],
            MEMBER["European Terrestrial Reference Frame 1993"],
            MEMBER["European Terrestrial Reference Frame 1994"],
            MEMBER["European Terrestrial Reference Frame 1996"],
            MEMBER["European Terrestrial Reference Frame 1997"],
            MEMBER["European Terrestrial Reference Frame 2000"],
            MEMBER["European Terrestrial Reference Frame 2005"],
            MEMBER["European Terrestrial Reference Frame 2014"],
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[0.1]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4258]],
    CONVERSION["UTM zone 32N",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Engineering survey, topographic mapping."],
        AREA["Europe between 6°E and 12°E: Austria; Belgium; Denmark - onshore and offshore; Germany - onshore and offshore; Norway including - onshore and offshore; Spain - offshore."],
        BBOX[38.76,6,84.33,12.01]],
    ID["EPSG",25832]]
st_crs(playgrounds)
Coordinate Reference System: NA
# Assign the CRS to playgrounds - its CRS is actually 4326, it just got lost in
# transport
st_crs(playgrounds) <- 4326

playgrounds <- playgrounds %>%
    st_set_crs(4326)

# Determine the CRS for the aarhus and elevation rasters
crs(aarhus)
[1] "PROJCRS[\"ETRS89 / UTM zone 32N\",\n    BASEGEOGCRS[\"ETRS89\",\n        DATUM[\"IRENET95\",\n            ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n                LENGTHUNIT[\"metre\",1,\n                    ID[\"EPSG\",9001]]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433,\n                ID[\"EPSG\",9122]]]],\n    CONVERSION[\"Transverse Mercator\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",9,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",0,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"easting\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]],\n        AXIS[\"northing\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1,\n                ID[\"EPSG\",9001]]]]"
crs(elevation)
Coordinate Reference System:
Deprecated Proj.4 representation: +proj=longlat +ellps=WGS84 +no_defs 
WKT2 2019 representation:
GEOGCRS["unknown",
    DATUM["Unknown based on WGS 84 ellipsoid",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1],
            ID["EPSG",7030]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8901]],
    CS[ellipsoidal,2],
        AXIS["longitude",east,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]],
        AXIS["latitude",north,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]] 
# Assign the CRS to elevation
crs_2 <- "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"
crs(elevation) <- 4326

Nice! You can determine what the CRS is using st_crs() for vectors or crs() for rasters. If one doesn’t exist, you can also use those functions to define the CRS.


Task 10: Transform your layers to a common CRS

In the previous exercise, when you ran st_crs() and crs() you may have noticed that the CRS’ were different for the different layers. The raster layer’s CRS began with +proj=longlat and the vector layer’s began with +proj=utm. In order to use these layers together in spatial analysis we will want them to have the same CRS.

In this exercise you will transform (sometimes this is called “project”) the objects so they share a single CRS. It is generally best to perform spatial analysis with layers that have a projected CRS (and some functions require this). To determine if your object has a projected CRS you can look at the first part of the result from st_crs() or crs() – if it begins with +proj=longlat then your CRS is unprojected.

Note that you will use method = "ngb" in your call to projectRaster() to prevent distortion in the elevation image.

Instructions

  • Use the crs() function with argument asText = TRUE on the aarhus layer to get the CRS as a string and save this as the_crs.
  • Use st_transform() to transform the vector playgrounds object to the CRS in the_crs.
  • Use projectRaster() to transform the raster elevation object to the CRS in the_crs. This will take a few seconds.
  • Use st_crs() on playgrounds and elevation to confirm that they now have the same CRS. They should all have a CRS that starts with +proj=utm.*

Solution

# Get the CRS from the aarhus object
the_crs <- crs(aarhus)
# , asText = TRUE)

# Project playgrounds to match the CRS of aarhus
playgrounds_crs <- st_transform(playgrounds, crs = the_crs)

# Project to match the CRS of aarhus
elevation_crs <- projectRaster(elevation, crs = the_crs, method = "ngb")

# Look at the CRS to see if they match
st_crs(playgrounds_crs)
Coordinate Reference System:
  User input: PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        DATUM["IRENET95",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]] 
  wkt:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        DATUM["IRENET95",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]
crs(elevation_crs)
Coordinate Reference System:
Deprecated Proj.4 representation:
 +proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs 
WKT2 2019 representation:
PROJCRS["ETRS89 / UTM zone 32N",
    BASEGEOGCRS["ETRS89",
        DATUM["IRENET95",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",9,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.9996,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",500000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Great work! This may be the least fun part of spatial analysis but knowing how to do it will save you from a lot of frustration later.

Task 11: Plot vector and raster together

If the layers do not share a common CRS they may not align on a plot. To illustrate, in this exercise, you will initially create a plot with the plot() function and try to add two layers that do not share the same CRS. You will then transform one layer’s CRS to match the other and you will plot this with both the plot() function and functions from the tmap package.

Note that for this exercise we returned all the layers to their original CRS and did not retain the changes you made in the last exercise.

With the plot() function you can plot multiple layers on the same map by calling plot() multiple times. You’ll need to add the argument add = TRUE to all calls to plot() after the first one and you need to run the code for all layers at once rather than line-by-line.

Instructions

  • Try to plot the playgrounds object on top of the aarhus object with plotRGB(aarhus) followed by plot(playgrounds, add = TRUE). Do you see the playgrounds? If not, use the st_transform to project to shared CRS.
  • Re-run the plot() code from the instruction above, color the playgrounds green and make the points stand out larger by setting the lwd argument to 3. You can also increase the transparency of aarhus, by setting its alpha at around 100.
  • Run the given tmap code. Note that tm_rgb() is used for multi-layered raster

Solution

# Plot aarhus and playgrounds (run both lines together) Do you see the
# playgrounds?
plotRGB(aarhus, stretch = "lin")  #, alpha = 100)
plot(playgrounds_crs, col = "green", lwd = 3, add = TRUE)

plot(playgrounds$geometry)

all.equal(st_crs(playgrounds_crs), crs(aarhus))
[1] "Modes: list, character"                              
[2] "names for target but not for current"                
[3] "Attributes: < Modes: list, NULL >"                   
[4] "Attributes: < names for target but not for current >"
[5] "Attributes: < current is not list-like >"            
[6] "Length mismatch: comparison on first 1 components"   
# Run the tmap code
library(tmap)
tm_shape(aarhus) + tm_rgb() + tm_shape(parks) + tm_polygons(col = "green") + tm_shape(playgrounds_crs) +
    tm_dots(col = "yellow", size = 1)

Great work! As you noticed, you mostly can’t plot layers together if they don’t have the same CRS. You’ll see later that there are exceptions but it is definitely best practice to ensure the layers you’ll work with have a common, projected CRS.

Task 12: The Detective

I have prepared a mystery.zip folder for you at this link: Download the mysterious files. They are all in Denmark, but they have different CRSs and file formats. Your job is to read them in R and create two plots: 1) one where you make all layers line up and 2) where you select and filter layers to create as pretty map as

If you need additional layers, check the Bonus section below

Bonus: Loading data from the internet

There are a lot of online resources for spatial data, such as OSM, AirBnB, GoogleMaps, DivaGIS, and GADM databases. Load some of this spatial data - vector and raster - using the geodata package. GADM has data for all current countries and is well-integrated with R. You can load the data directly into R with the function geodata::gadm() and start processing. The available data are:

  • SRTM 90 (elevation data with 30s - 3s resolution between latitude -60 and 60)
  • World Climate Data (Tmin, Tmax, Precip, BioClim)
  • Global adm. boundaries (different levels)

In the case of GADM you must also provide the level of administrative subdivision (0=country, 1=first level subdivision)

Instructions

  • To read elevation with the elevation_30s() function from the geodata package, select the following arguments:
    • country needs a 3 letter ISO code for Denmark; getData(‘ISO3’) let’s you see the codes, you can also specify Lat and Long arguments instead of country, if you need a specific area that overlaps multiple countries.
    • mask needs to be set to TRUE as we wish to mask surrounding countries.
  • to read administrative data with gadm()function, choose the following:
    • country needs a 3 letter ISO code for Denmark; getData(‘ISO3’) let’s you see the codes,
    • level should be set to the level of administrative subdivision (0=country border, 1=first level subdivision, eg. Midtjylland region, 2 = municipalities,… ).

Examples

library(geodata)
`?`(gadm())
# elevation <- elevation_30s(country = 'DNK', mask = TRUE, path = tempdir())
# municipalities <- gadm(country = 'DNK', level = 2, path = '.') %>% st_as_sf()
# DK_border <- gadm(country = 'DNK', level = 0, path = '.')%>% st_as_sf()